Dynamical current-current correlation of the hexagonal lattice and graphene 
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We discuss the dynamical current-current correlation function of the hexagonal lattice using a 
local current operator defined on a continuum-replica model of the original lattice model. In the 
Dirac approximation, the correlation function can be decomposed into a parallel and perpendicular 
contribution. We show that this is not possible for the hexagonal lattice even in the Dirac regime. 
A comparison between the analytical isotropic solution and the numerical results for the honeycomb 
lattice is given. 
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I. INTRODUCTION 

Graphene is a two-dimensional carbon allotrope which 
was isolated in 2004 1 and has attracted immense re- 
search activities due to its novel mechanical and elec- 
tronic properties. 2-5 Whereas the mechanical properties 
are determined by electrons with sp 2 -hybridization, the 
electronic properties can be mainly deduced considering 
only the 7r-electrons. The simplest model to study the 
electronic response of graphene to an external field or po- 
tential is thus given by a one-orbital tight-binding model 
on a hexagonal lattice. 

Most of the novel electronic properties of graphene 
originate from the fact that there are two equivalent 
atoms in the Wigner-Seitz cell which give rise to two 
gapless bands with linear density of states close to the 
neutrality point. Most standard results of solid state text 
books can thus not be applied to the case of graphene due 
to the different dispersion and/or dimensionality, but also 
due to the two coherently coupled bands. 

An example is the density-density correlation or Lind- 
hard function which in the case of the honeycomb lattice 
is given by 6 
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with the eigenenergies E ± (k) — ±t\<fi(k)\ (t « 2.7eV is 
the hopping amplitude), np(E) the Fermi function, g s — 
2 the spin-degeneracy and </>(k) the complex structure 
factor defined below. Due to the two gapless bands, the 
above expression contains the band-overlap function 
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which marks the crucial difference to the standard text- 
book results containing only one band. 7 

In the linear (Dirac) approximation of the band dis- 
persion, the above expression can be solved analytically 



for finite chemical potential \i at zero temperature. 8,9 In 
the static case, it shows differences to the one-band result 
by Stern 10 for |q| > 2kp with kp the Fermi wave vector 
due to the contribution of interband processes. For finite 
frequencies, these differences are even more pronounced 
and lead to a logarithmic singularity at hu = 2/i. 

The density-density correlation function or polariz- 
ability of graphene was calculated in a number of pa- 
pers using different formalisms and introducing various 
modifications to the original Dirac Hamiltonian. 11-23 Us- 
ing these results, plasmons, 8 > 9 > 20 > 24 wrinkles, 25 van-der- 
Waals interactions 26 and forces due to moving external 
charges 2 were discussed. In this paper, we will focus on 
the related current-current correlation function 7T 1 J (q, uj) 
of graphene = x,y) starting from the tight-binding 
model of the honeycomb lattice. 

In the Dirac approximation, the system is rotationally 
invariant and the current-current correlation function can 
be decomposed in a parallel 7r" and perpendicular con- 
tribution 7T . The parallel contribution is related to the 
density-density correlation function via the continuity 
equation and thus determines the dielectric properties of 
the system. The perpendicular contribution is related to 
the magnetic susceptibility which in the static case has 
been first discussed by McClure 28 via the Hclmholtz free 
energy and recently by Ando and co-workers using n , 29 
In view of new experiments on the magnetic behavior of 
graphene 30 , the magnetic susceptibility was also calcu- 
lated including electron-electron interactions to first or- 
der which results in a paramagnetic response away from 
half-filling. 31 

Here, we shall mainly discuss 7r J -(q, uj) for finite fre- 
quencies. In the Dirac approximation, this was first 
done in Ref. 32 . We will summarize their results and 
compare the analytical solution of the isotropic system 
with the numerical solution of the hexagonal lattice. For 
that, we will define a local current operator defined for a 
continuous-replica model of the original lattice Hamilto- 
nian. This formalism permits deeper insight in the lattice 
effects and can be used to calculate corrections which are 
lost in the scaling limit, i.e., the Dirac model. 

The paper is organized as follows. In section II, we will 



2 



define the continuum model and derive the local current 
operator of this model. We will further show that this op- 
erator satisfies the continuity equation with respect to the 
density operator defined on the lattice. In section III, we 
will present general expressions for the current-current 
correlation function and introduce the parallel and per- 
pendicular contribution defined for the Dirac model. In 
section IV, we summarize the analytical results and com- 
pare them with the numerical results obtained from the 
hexagonal lattice. We close with a summary and con- 
clusions and give real expressions for the current-current 
correlation function in an appendix. 



II. CONTINUUM MODEL AND CURRENT 
OPERATOR 

To calculate the current-current correlation function 
for a lattice model for finite wave vector q, we are con- 
fronted with the following problem. The current operator 
for a lattice model, as given by the continuity equation, 
describes the flow from site i to j per unit time. 33 In or- 
der to define a vector which depends on one lattice site 
instead of two, one needs to define a continuous model 
based on the Hamiltonian in reciprocal space. If the vec- 
tor potential A is a smooth function of r, then the cou- 
pling between A and the current can only see the smooth 
part and the continuous limit is justified. 

Let us start with the tight-binding Hamiltonian of 
a general bipartite lattice with N c lattice sites R and 
nearest- neighbor lattice vectors 5: 



R,<5 



(3) 



The spin-index on the operators shall be suppressed 
throughout this work. With the Fourier components 
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this reads for £r.r +< 5 = t 
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(4) 
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with 0(k) = J2 6 



ik<5 



the complex structure factor where 



the sum goes over all nearest-neighbor vectors 5. Notice 
that the phase factor e lk<5 in Eq. (5) is important for the 
definition of the current. 34 

We will now define a continuous model by introducing 
the following Fourier components 
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where c = a,b and A the area of the sample. 

The continuous version of the Hamiltonian thus reads 



H 



= -t j d 2 r [at 



(r)<£(-iV)b(r) + H.c 



(9) 



The gauged Hamiltonian is obtained by replacing — iV — > 
— iV + f A(r) (e > 0). Notice that by going back in 
Fourier space, we obtain the correct Peierls substitution 



t R ,R+s -> t R M + se^ & +t dlA « 



(10) 



in the case of a gauge field which is constant over one 
lattice spacing, i.e., for a spatially weakly varying field. 
Because e v s b(r) = b(r + S), we can write Eq. (9) as 

H=-t f d 2 rJ2 K(r)fo(r + S) + H.c] . (11) 

The continuous model thus consists of infinitely many 
replica of the original lattice model. The unperturbed 
Hamiltonian is homogeneous, real (not crystalline) mo- 
mentum is conserved and yet, each particle is bound to 
hop in the replica where it lives with strict fidelity to the 
original lattice Hamiltonian. In particular, the lattice 
anisotropy is fully preserved. Also, the minimal substi- 
tution used to include the perturbing vector potential 
guarantees gauge invariance to all orders. 

For this model, the current can be defined by 



j(r) 



SH 
~6A(r) 



= j P (r)+j D (r)+0(A 2 ), (12) 



where the diamagnetic contribution j D is linear in the 
gauge field A. 

For the paramagnetic operator, we obtain 

j p (r) = ^ V S f ds [at(r - sS)b(r + (1 - s)S)] + H.c. 

(13) 

which consists of a symmetrized version of the paramag- 
netic current given in Refs. 35,36 which is obtained from 
the above formula by setting s = 0. 

For the diamagnetic contribution, we obtain with 
(summation over j is implied) 



' 4 (r) = J d 2 r 
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A=0 



the following expression: 



J D ' 4 ( r ) =~ t iyi ^ / dsds ' t at ( r - sS ) b ( r + (! - s )<*) 
h x Jo 



x {sA j (r - ss'S) + s'A j (r + ss'S)} + H.c] 



(15) 



which again resembles a symmetrized version of the dia- 
magnetic current given in Refs. 35,36 which this time is 
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obtained from the above formula by setting s = and 
s' = 1. Notice that the diamagnetic current is non-local 
in the external gauge field. 

In linear response, only ground-state averages enter in 
the diamagnetic current. With the energy per bond per 
unit area 



/ibo„d - -2t(a^(r) b(r + S)) 



(16) 



which is independent of both r and S, the Fourier trans- 
form of the paramagnetic and diamagnetic current are 
given by: 

tf* ^Ef(k,q)«K+q + (^k-q)) , ^k +q , (17) 



h 

with 



ik+q)-<5 _ 



(18) 

(19) 
(20) 



We obtain the general expression for the current- 
current correlation function 
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(22) 



with E ± (k) = ±i|^(k)| and n F (E) = (e^ E ~^ + l)" 1 
the Fermi function. 

This is the same expression as for the density-density 
correlation function of Eq. (1), but the band-overlap is 
now given by 



(23) 
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For q — > 0, we obtain the same expression as in Refs. 36,37 . 
The Fourier transform of the particle density of the 
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lattice model is given by n q = X)k( a k a k-K 
For the charge density p q = en q , the continuity equation 
Pq — iq • j q — is obeyed for the paramagnetic current 
operator of Eq. (17). We can thus consider this operator 
to be the current operator of the lattice model for general 
q. In the same manner, the diamagnetic term is also 
correct for arbitrary q. 



Due to charge conservation, we have 
e 2 u; 2 Im7r ' (q, to) — qiqjlmir 1 ' 3 (q, u>) where summa- 
tion over double indices is implied. To see this within 
our notation, we note that gi</>*(k, q) = (f>(k) — <fi(k + q) 
and thus 



/£ (k,q)= g V 



/± J (k,q) 



(|^(k)| T |0(k + q)|) 2 



(24) 



III. CORRELATION FUNCTION 

We can now determine the current-current correlation 
function. In terms of the bosonic Matsubara frequencies 
hw n = 2nn//3 (/3 = l/k B T), it is defined by 

1 r h ' 3 

^(q.iwn) = j^J o dre^^ir)^) . (21) 



which proves the relation since (ftu/t) 2 — (|</>(k)|=F |</>(k+ 
q)l) 2 - 

In the Dirac cone approximation, the expressions sim- 
plify considerably. Denoting the angle between k and q 
by p and neglecting terms proportional to sintp which 
cancel to zero due to the angle integration, we have for 
the effective band overlap 



1 /3a 
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q 2 |k||k + q; 



q 2 |k||k + q| 
I 



2 Q 
1 — 2 sin p + — cos 03 
k 



1 — 2 sin p + — cos ip 



k 



(25) 
(26) 



where we introduced the carbon-carbon distance a = 
0.14nm. 



The system linearized around the Dirac point is ro- 
tationally invariant. We can thus decompose 7r lJ into 
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a longitudinal component 7r" and transverse component 
n- 1 . These are defined by Eq. (22) after substitution of 
the overlap function f± by 



IK_L) _ 1 /3a \ / k + qcosip - 2£;sin (p 

f± "2 2 l 1±(T) jk+^j 



where the energy per bond per unit area of the hexagonal 
lattice is given by 



h ho n d = |jE S+ ( k ) [n F (E-(k))-n F (E+(k))} 



(27) and the band cutoff = 6t. 



(34) 



We then recover the general relation 

^'(q.u;) = S*"(|qM + (*J SK^qM 



q- 



We note that the overlap function f±° in the Dirac ap- 
proximation is proportional to /J. , but with the last term, 
2k sin 2 ip, missing. 8 

Due to current conservation and q 2 n^ — qiir^qj, the 
parallel component of the current-current correlation is 
related to the density-density correlation by 

q 2 J (\q\,w) = -<[p q , q • j_ q ])/(ftA) + e 2 u,V>°(|q|, uj) . 

(29) 

Apart from the constant surface or contact term, which 
was determined in Ref. 39 for the linearized Dirac model, 
we are thus left with the calculation of the perpendicular 
component ir 1 - which is related to the magnetic suscepti- 
bility XM(q,w)//i = ir^(c\,u)/\q\ 2 for uj < |q| with ^ 
the magnetic permeability. 33 

For the full dispersion, we have — ([p q , q- j_ q ]) /(hA) = 
QiXq ' 1 '^ Qj- K is thus often more transparent to deal with 
the physical response, IT'- 7 , which includes the diamag- 
netic contribution: 



(28) 



ir^q, W )=7r'''(q,w)+X?^ 



(30) 



Charge conservation then implies 

q i Il^(c l ,uj)q j = e 2 uj\ '°(q,uj) . (31) 

Notice that the anisotropy of the response for finite q 
requires the full tensorial structure of IP'- 7 . In particu- 
lar, the relation between polarizability and conductivity 
reads 



a; a % ' 3 {q,uj) qj = iuje 2 Tr 0,0 (q, uj) . 



(32) 



We will show in the next section that the often used 
scalar version of Eq. (32) would not hold for the lattice 
model even in the regime where the Dirac approximation 
is justified. 

We finally state the general f-sum rule for a bipartite 
tight-binding model: 

dujujlmn ' (q,uj) = l^^sin 2 (^) (33) 



IV. RESULTS 

We will now summarize the analytical results obtained 
for the Dirac approximation at zero temperature first pre- 
sented in Ref. 32 and compare them with the numerical 
results obtained from the hexagonal lattice. 



A. Analytical results 

In order to present the analytical results, we express 
the current-current correlation function of Eq. (22), 
7r ± (g, uj), by two dimensionless functions 



eH 

h 2 



[*£(q,u>) + Arf(q,wj\ , (35) 



where we will use the the superindex + to denote the lon- 
gitudinal component (||) and the superindex — to denote 
the transverse component (_L). We restrict the discus- 
sion to uj > since n ± (q, —uj) = [n ± (q,uj)] and to fi > 
due to particle-hole symmetry. 7r^ contains the contri- 
bution for the system at half-filling, i.e., interband con- 
tributions, whereas Air^ contains the contributions due 
to the finite chemical potential /i, i.e., intraband contri- 
butions. The formulas are given in terms of the Fermi 
velocity hvp — |ai. 

The results can be written in compact form using two 
dimensionless, complex functions defined as 



ti+/ \ q fujj 



i 



G ± (x) = xVx 2 — 1 =F In (x + vV - l) 



(36) 
(37) 



Let us first present the results for the undoped system. 
For large energy cutoff A^ ^> 1, we have 



_9^Ae 
8tt t 



mF ± (q,uj) 



(38) 



Notice that the constant cutoff term can be obtained ei- 
ther from the Kramers-Kronig relation or from the conti- 
nuity equation. This connection gives rise to the so-called 
f-sum rule. 39 

The contribution due to the finite chemical potential 
reads 



An±(q,uj) = ±^ 



- 2nt{vFq) ^^){GH X+ ) (39) 
-6 (jc_ - 1) [G ± (x_) T iTr] - 6 (1 - x-) G ± {-x-)} 



where we defined x± 



hv F q 

The above expression for graphene shall be contrasted 
with the expression for the two-dimensional electron gas. 
For quadratic dispersion ek = S 2 k 2 /(2m), we have 



lqla-0. 1, ^/t=0.05 



'few) 




2tt q 2 2?r 



21 T 2 




where the term proportional to /i = £fc F corresponds to 
the contact term which is canceled by the diamagnetic 
contribution. 

Eq. (39) can be written as real and imaginary part in 
terms of three real dimensionless functions 



g huj 



1 



(41) 



G>(x) = xs/ x 2 — 1 =F cos h 1 ( x ), x>l, 
1 — cos _1 (a;) , |x| < 1 . 



G%{x) = ±xVT 



The lengthy expressions are given in the appendix. 

Let us now discuss two limiting cases. For the long 
wavelength limit q — > 0, we obtain 



2/i 


V In 


2^ + 






2^i - 





i-e(hw - 2,u) 



(42) 



Using the RPA-approximation for the longitudinal part, 
the above expansion leads to plasmon excitations for 
which the logarithmic term is usually neglected. 8 ' 9 Due 
to the sign change of the photon propagator in the 
case of transverse modes, the denominator of the RPA- 
approximation cannot become zero for the perpendicular 
part without the logarithmic term. But including it leads 
to a new transverse electromagnetic mode in graphene. 40 
For the static case, we obtain the following formula 
which was first given in Ref. 29 : 



7T (q,u = 0) 



9 -v F qe(q-2k F )G~(^) (43) 

q 



h 8tt 



where kp = /i/ (hvp). The parallel component ir + is zero. 
For fixed q, tt~ is only non-zero for /i < Ympqjl and since 

j} dxG<(x) = 4/3, the limit q — >• leads to the well 
known delta function for the diamagnetic susceptibility 
of graphene: 



XM = -(J,a — e 2 v F 5(n) 

D7T 



Numerical Results 



(44) 



We shall now compare the analytical results of the lin- 
earized, isotropic Dirac model with the numerical results 
obtained from the hexagonal lattice. In Fig. 1, we show 
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FIG. 1: The imaginary part of the current-current correla- 
tion function Ini7r*' l (qa;, q y , w) as function of the energy uj at 
k B T/t = 0.01 for different directions with |q|a = 0.1. The 
results obtained from the Dirac-cone approximation 7r"(g, w) 
and ir ± (q,uj) are also shown (dashed lines). 



the imaginary part of the current-current correlation 
function Iunr' t '' l '(q x , q v , uS) as function of the energy fko at 
ksT/t = 0.01 for different directions with |q|a = 0.1. 
The results obtained from the Dirac-cone approxima- 
tion n"(q, u)) and TT ± (q,u}) are also shown (dashed lines). 
Clearly, there are strong differences for energies hu > t 
due to the van Hove singularity. The inset shows that 
there is a peak splitting for the different directions due 
to the different contributions of the three M-points, also 
present in the charge response. 21 

In Fig. 2, the same curves are shown, but for lower en- 
ergies. On the left hand side, the wave vector q is parallel 
to the current and on the right hand side perpendicular. 
The results obtained from the Dirac-cone approximation 
w"(q,u) and TT ± (q,uj) are also shown (dashed lines). 

For the perpendicular contribution of 7r M (right hand 
side), clear differences are seen for lower energies due to 
the finite temperature ksT/t = 0.01 used in the numer- 
ical calculation. This results in a thermal broadening of 
the delta-function of Eq. (44) and is responsible for the 
diamagnetism found in graphene 30 since intrinsic dop- 
ing leads to [i ^ 0. The dotted lines on the right hand 
side refer to the same curves but at lower temperature 
kBT/t = 0.001 which agrees well with the Dirac cone 
approximation now also at low energies. 

When q is in ^-direction which was chosen to be the 
high symmetry axis which connects the T- and M-point 
of the Brillouin zone, there is good agreement with the 
result coming from the Dirac cone approximation (except 
for the deviations in tt 1 - due to temperature, mentioned 
before). When q is in ^/-direction, we observe a peak 
splitting around the resonant energy hu: = hvpq (see in- 
set on the left hand side). There are thus lattice effects 
which show up even in the regime where the Dirac cone 
approximation and where the system should be isotropi- 
cally invariant. 
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lqla=0.J , jl=0.05t 



\q\a=0.1 , n=0.05t 




FIG. 2: The imaginary part of the current-current correla- 
tion function Ini7r 1,l (g a; , q y , oj) as function of the energy uj 
at ksT/t = 0.01 for different directions with |q|a = 0.1. 
Left hand side: The wave vector q is parallel to the cur- 
rent. Right hand side: The wave vector q is perpendicular to 
the current. The results obtained from the Dirac-cone approx- 
imation 7r"(q,u;) and ir ± (q,uj) are also shown, respectively 
(dashed lines). Inset: Energy region around the resonant en- 
ergy uj — VFq- Right hand side: Also the curves for lower 
temperature ksT/t — 0.001 are shown (dotted lines). 



V. SUMMARY AND CONCLUSIONS 

We have discussed the dynamical current-current cor- 
relation function of the hexagonal lattice and of graphene 
modeled by the linearized Dirac model. To define a local 
current operator, we introduced a continuum-replica of 
the original lattice model. The resulting paramagnetic 
current operator obeys the continuity equation with re- 
spect to the density operator defined on the original lat- 
tice. The diamagnetic response is non-local. 

We then gave explicit expressions of the current- 
current correlation function for the honeycomb lattice 
and defined the longitudinal and transverse component 
in case of the rotationally invariant Dirac model. For the 
Dirac model, explicit analytical expressions were given 
where the results for the longitudinal component can be 
obtained via the continuity equation from the density- 
density correlation function, as was discussed in detail. 

In the last part of this paper, we showed that in the 
honeycomb lattice, the longitudinal and transverse com- 
ponent cannot be defined for energies around the res- 



onant energy huj — hvpq. This is reminiscent to the 
fact that also the polarizability is not well described by 
the Dirac approximation for these energies. 21 The scalar 
relation between the conductivity and the polarizability 
which makes use of the fact that there is a parallel compo- 
nent does thus not hold for the lattice model. This might 
be important for first principle studies which make use 
of this relation. 
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FIG. 3: Display of the different regions characterizing the 
current-current correlation function given in the appendix. 
The regions are limited by straight lines ui = q (solid), ui = 
q — 2/i (dashed) and u = 2fi — q (dotted) where we set h = 
vf = 1. 
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VII. APPENDIX 



Here, we shall present the real and imaginary part of 
7T^~ in terms of the three real dimensionless functions 



J KH ' ' 16tt t 



G>(x) — x\/ x 2 — 1 =F cosh 1 (x), x > 1 , 
G%{x) = ±xyjl - x 2 - cos^x) , |x|<l. 



For the imaginary part, the additional terms at finite 
doping then read in the language of Fig. 3: 



ImA7T^(g,w) 



>^ hvpq ' >\ tarn 



-f (q,u) x < 



£j± / 2fi-\-huj \ 



-G<( 


I o 



hvpq 



1 A 

1 B 

2 A 

2 B 

3 A 
3 B 



(46) 
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For the real part, we get in the language of Fig. 3: 



2nt (vpq)'- 



hv F q I 

s-i± ( hm — 2fi \ j_ ynf± ( 2[i+huj \ 

Hv F q } f <^ hv F q ) 

>\ hvpq ' >^ hvpq ' 



1 A 

1 B 

2 A 

2 B 

3 A 
3 B 



(47) 



Since this agrees with the complex expression given in Eq. (39) . 

For more details, see Ref. 8 and Ref. 32 . 



G± w={±il> :w<i (48) 
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